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We report on Monte Carlo studies of the kinetic exchange 
model for (III,Mn)V ferromagnetic semiconductors in which 
5 = 5/2 local moments, representing Mn 2+ ions, are exchange 
coupled to band electrons. We treat the Mn 2+ spin orien- 
tations as classical degrees of freedom and use the Hybrid 
Monte Carlo algorithm to explore thermodynamically impor- 
tant Mn spin configurations. The critical temperature T c of 
the model is unambiguously signalled in our finite-size simu- 
lations by pronounced peaks in fluctuations of both Mn and 
band carrier total spins. The T c 's we obtain are, over much of 
the model's parameter space, substantially smaller than those 
estimated using mean-field theory. When mean-field theory 
fails, short-range magnetic order and finite local carrier spin 
polarisation are present for temperatures substantially larger 
than Tc. For the simplest version of the model, which has a 
single parabolic band with effective mass m* , the dependence 
of T c on m* is sublinear at large masses, in disagreement with 
the mean-field theory result T c tx to* . 

PACS numbers: 75.70.Dd, 75.40.Mg, 75.40.Cx 



I. INTRODUCTION 

Because their electronic properties are, at low temper- 
atures, extremely sensitive to external magnetic fields, 
diluted magnetic semiconductors have long attracted at- 
tention Q]. The discovery of ferromagnetism at com- 
paratively high temperatures (T c > 100K) in Mn- 
doped (III,V) semiconductors Hi| has intensified inter- 
est, partly because of the new technological pathways 
that might be opened if room temperature ferromag- 
netism were achieved in a semiconductor with favorable 
materials properties The search for materials in 

this class with higher ferromagnetic transition tempera- 
tures, and the effort to achieve a deeper understanding of 
all physical property trends, has stimulated a large body 
of experimental and theoretical research. 

To date, most theoretical work |^-||] has been based 
on kinetic exchange models of the type used here and de- 
tailed below. First principles density- functional-theory 
electronic structure calculations |ic],[ll| are qualitatively 
in accord with these phenomenological models, although 
more work may be necessary to reliably judge the ac- 
curacy of the local-spin-density approximation for these 
materials before they can be used to guide model re- 
finement. The model's low-energy degrees of freedom 
are S — 5/2 spins representing Mn 2+ ions and holes in 



the semiconductor valence band that are free to move 
through the system. Ferromagnetism is most easily dis- 
cussed J|-§1 in a Mn continuum version of the model, jus- 
tified in part by the small ratio of hole density to Mn ion 
density in experimental systems. The constant Mn den- 
sity model is somewhat analogous to the jellium model 
for metals in that it artificially preserves translational 
invariance. In modelling these materials, the simplest 
plausible approximation is one that both represents lo- 
cal moments by a uniform density continuum and ignores 
correlations between local moment orientations and the 
state of the valence-band-hole gas. In the following we re- 
fer to this combination of idealizations as the mean-field- 
approximation. In the mean-field-approximation 
the ground state is always ferromagnetic because the ex- 
change energy gained by polarizing is linear in the hole- 
spin density whereas the band and hole-hole interaction 
energy cost is quadratic. For the Mn concentrations at 
which the highest ferromagnetic transition temperatures 
occur, it appears that mean-field-theory predictions for 
T c , magnetic anisotropy, optical absorption, and other 
properties are at least qualitatively in accord with ex- 
periment, although further experimental and theoretical 
work is needed to confidently establish where this ap- 
proach is and is not reliable. 

It is well known that mean-field theories tend to overes- 
timate the stability of ordered phases, because they fail 
to account for thermal and quantum fluctuations that 
can destroy long-range order. In particular, they pre- 
dict critical temperatures which are systematically too 
high. A step toward a more complete description of 
ferromagnetism in diluted magnetic semiconductors was 
taken in Ref. |0] , where a theory of collective spin- wave 
excitations was developed for the continuum version of 
the kinetic-exchange model. T c estimates jOl based on 
this theory of the ferromagnet's spin stiffness establish 
that mean-field theory becomes less reliable for stronger 
exchange coupling and for flatter semiconductor valence 
bands; the higher the predicted T c the lower its relia- 
bility. The approximate collective-fluctuation T c bounds 
proposed in Ref. |l3| are based on T — spin-wave en- 
ergies of the continuum version of the kinetic exchange 
model. The considerations of this paper demonstrate 
that collective magnetization fluctuations will control the 
temperature at which long range magnetic order is lost 
in at least a part of the model-parameter range relevant 
to ferromagnetic semiconductors. In this work we report 
on Monte Carlo simulations, based on a Hybrid Monte 
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Carlo algorithm ]TJJ| , that fully account for thermal fluc- 
tuations of Mn spin orientations and shed further light on 
the finite-temperature properties of the kinetic-exchange 
model. Some of the results of these calculations were 
described briefly in Ref. |l3|]. In this paper we give a 
complete description of the calculations, including all im- 
portant technical details. In addition, we also present 
results obtained using a realistic six-band k ■ p model for 
the valence bands of zincblende semiconductors, allow- 
ing us to assess the importance for thermal fluctuations 
of band structure details. A related Monte Carlo study 
by Sakai et al. fl5|| has the same motivation as ours, but 
examines a model of Ising localized spins coupled to band 
electrons. This model has gapped collective excitations 
and qualitatively different thermodynamic properties. 

Our simulations are directed toward the regime of 
Mn concentration in IIIi_ a; Mn a ;V alloys where the high- 
est ferromagnetic transition temperatures have been ob- 
served, x ~ 0.05 |^]. In this regime the heavily- 
doped semiconductors are strongly disordered three- 
dimensional metals with kp£ ~ 3. (Here Icf is a typical 
Fermi wavevector and £ is a mean-free path estimated 
from the measured T = resistivity.) The sources of 
disorder in these alloys have not been fully cataloged, 
although randomness in the Mn alloy positions and an- 
tisite defects (group V elements on group III sites) are 
certainly among them. In our simulations we include 
only the former disorder source. In the continuum ver- 
sion of the model this source is absent and the system 
is disorder-free. We find that disorder has quantitative 
but not qualitative importance for the model we study. 
Randomness in the Mn ion distribution will be increas- 
ingly important at lower Mn concentrations, particularly 
on the insulating side of the metal insulator transition 
which occurs for x ~ 0.02, and has an overriding impor- 
tance in the models studied in Ref. Jl6| that are directed 
to the dilute Mn concentration limit. 

In the extreme dilute limit, all holes are bound to Mn 2+ 
acceptors, forming neutral complexes. For small Mn con- 
centrations, holes may hop from acceptor site to acceptor 
site, but it is still useful to regard them as belonging to an 
impurity band. On the metallic side of the metal insula- 
tor transition, the concept of an impurity band is less use- 
ful and the host semiconductor band structure has to be 
modelled realistically, as in the present work. In order to 
correctly capture the physics of the dilute limit, we would 
need to include Coulombic interactions between Mn ions 
and valence band spins in addition to exchange interac- 
tions. If hole-hole interactions were treated in a mean- 
field-theory, including Coulomb scattering would add a 
self-consistently determined scalar potential to the band- 
electron Hamiltonian. There are several obstacles to in- 
cluding this scalar potential in our work at present. The 
most fundamental problem is current ignorance about 
the distribution of charged defects in these materials. 
It is known that the hole density in typical samples is 
approximately three times smaller than the Mn density 



and suspected that the compensation is due mainly to 
the antisite defects present in low-temperature-growth 
compound semiconductors |l7| ]. These defects will play 
as important a role as the Mn ions in determining the 
screened scalar potential. There are also technical ob- 
stacles to including Coulomb interactions in our Monte 
Carlo calculations. If they were explicitly included, we 
would have to self-consistently solve for the valence band 
electron screening charges for each Mn spin configura- 
tion, or adopt a Car-Parrinello [Q type of scheme. In 
the present study we shall therefore ignore Coulomb in- 
teractions, partially to enable progress on a reasonable 
level of technical difficulty warranted by the present sta- 
tus of experimental insight and microscopic theoretical 
modelling. This approximation dominates so far the in- 
terpretation of experimental results on magnetic proper- 
ties of Mn-doped III-V semiconductors, and as well the 
theoretical discussion. 

This paper is organized as follows: In section || we 
discuss the kinetic-exchange model of Mn-doped (III,V) 
semiconductors. In section [II we summarize the main 



features of the Hybrid Monte Carlo algorithm that are 
important for our purposes. Numerical results obtained 
by this method as well as a detailed discussion of im- 
portant technical aspects are presented in section IV. 



Our data show very clearly the typical features of ferro- 
magnetic transitions signalled unambiguously by suscep- 
tibility peaks. We discuss the critical temperature as a 
function of effective mass, carrier density, and exchange- 
coupling strength. We close with a summary in section 
V. 



II. THE MODEL 

We now describe the kinetic-exchange model we use 
for (III,Mn)V diluted magnetic semiconductor ferromag- 
nets in detail. In this material class, the highest critical 
temperature of 110K has been obtained for Gai-^Mn^As 
with x fa 0.05 [§. The manganese ions are assumed to 
be Mn 2+ with a 3<i 5 , S = 5/2 configuration, and to be 
placed at random on group III clement sites so that they 
are acceptors. The concentration of holes is, however, 
allowed to differ from the Mn concentration, since the 
system is strongly compensated by anti-site defects, i.e. 
group V atoms on group III sites. In practice the hole 
concentration is often not accurately known; its determi- 
nation from Hall effect measurements is complicated by 
the extraordinary Hall effect that occurs in all ferromag- 
netic metals 

Long range magnetic order occurs in this dilute system 
of magnetic moments because of an exchange interaction 
between the spins of the manganese ions and spins of the 
valence band holes. The minimal model to describe this 
carrier-induced ferromagnetism is 



n = E £ + E / d * rJ ^- RxW) ■ sj 



(i) 
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This Hamiltonian describes noninteracting carriers in a 
parabolic band characterized by an effective mass m*, 
whose spin density s(f) is coupled to manganese ion spins 
Si at locations Ri by an antiferromagnetic exchange cou- 
pling. To account for the finite extent of the manganese 
ions, the exchange is modelled by a spatially extended 
coupling J(r) |19fl ; we use a Gaussian centered around 
the ion position, 



J(r) 



pd 



(2) 



Both the strength and range of this interaction are phe- 
nomenological parameters to be fixed by comparison with 
experiment or, ideally, to be extracted from first prin- 
ciples electronic structure calculations. The exchange- 
coupling range parameter ao in Eq. (^) is required in our 
calculations to keep exchange-coupling shifts in quasi- 
particle energies finite. In the numerical studies to be 
presented below we usually choose ao = O.lnm, which 
appears to be a reasonable value for the d-ion radius in 
Gai_ x Mn x As, whose lattice constant is ~ 0.60nm p||3||. 
However, we will see that properties of the system de- 
pend only very weakly on ao, as long as this quantity is 
not increased to large values. 

The exchange parameter J p d > characterizes the 
strength of the exchange interaction between carriers and 
manganese ions. Typical experimental values for this 
quantity are of the order of O.leVnm 3 . (The pd subscript 
on the coupling constant is a reminder of its presumed 
physical origin in the hybridization of valence-band p- 
orbitals and Mn-ion d-orbitals. Since the carrier spin 
density s(r) has dimension 1 /volume and Si is dimen- 
sionless, J p d has the dimension of energy times volume.) 

We emphasize that the interaction in Eq. ([!]) is of the 
isotropic Heisenberg type. In a parabolic band model, 
total energies are therefore invariant under global rota- 
tions of all Mn spins and slow variations in orientation 
have a small energy cost. Even with the realistic band 
models discussed below, the energy cost of global Mn spin 
rotations is still very small . Models with Ising-like in- 
teractions |H| do not capture this essential element of 
(III,Mn)V ferromagnet physics, although they are much 
more easily studied. 

An important refinement of the above model, Eq. (Q), 
is obtained by replacing the single parabolic band by a 
realistic model for the host semiconductor valence bands. 
In the materials of interest these derive from atomic p- 
orbitals and have a spin-orbit coupling strength that ex- 
ceeds the Fermi energy. Importantly, eigenstates do not 
have definite angular momentum quantum numbers at 
finite wavevector; at the center of the Brillouin zone, 
however, they have definite total angular momentum 
J G {3/2, 1/2}. We use these states in a k ■ p scheme, 
choosing as basis states p0|,|7[| 
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Near the Brillouin-zone center, the band structure can 
be parameterized by a small number of symmetry- 
adapted parameters of the Kohn-Luttinger Hamiltonian 
fjj"| , po| j7| ) . In the above basis this effective kinetic-energy 
operator reads 
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where the entries of this 6x6 matrix for each wavevector 
k are given by 



Hhh = ^ [(7i + I2)(kl + k 2 v ) + (71 - 2 72 )fc z 2 ] 

Hlh = i [ (71 " 72)(fc ' + k y ] + (71 + 272)fc '] 



2m 



m 

2m 
V2h 2 
2m 



[l2(k 2 x - k^) - 2i l3 k x k y ] 
l2 [2k 2 z -(k 2 x + k 2 y )] . 



(5) 



IU) 



The quantities 71, 72, 73 are called Luttinger parame- 
ters, m is the bare electron mass, and A so is the spin- 
orbit coupling which splits the six states at the valence 
band edge into a doublet and a quartet. These band 
parameters are accurately known for a large number of 
(III,V) compound semiconductors P^| . For all numeri- 
cal calculations we use (71,72,73) = (6.85,2.1,2.9) and 
A so = 0.34eV, the band parameter set appropriate for 
GaAs. 

In the kinetic-exchange model with realistic band 
structure, the Hamiltonian (^) replaces the kinetic part of 
Eq. (|). The Hilbert space dimension of the one-particle 
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problem is, thus, enlarged by a factor of three. The 
Kohn-Luttinger Hamiltonian has already been adopted 
in some mean-field calculations for the continuum ki- 
netic exchange model H|t]|. This realistic description 
is essential if important magnetic parameters connected 
with magnetocrystalline anisotropy and anisotropic mag- 
netoresistance are to be addressed. 

The material Gai-^Mn^As has a zincblende structure, 
in which the cations and anions each form an fee lattice. 
In the alloy manganese ions are located on randomly 
chosen cation sites, forming a disordered system of lo- 
cal moments. In an approximation consistent with our 
band Hamiltonians we ignore the underlying lattice and 
place the ions positions completely at random within pe- 
riodically continued cubic simulation cells. In the Monte 
Carlo calculations described below, we average over dif- 
ferent disorder realisations. However, we find that the re- 
sults of different realizations typically lie within a range 
of a few percent. 

In our simulations we treat the S = 5/2 Mn spins as 
classical degrees of freedom. This approximation is not 
severe even in a mean-field theory where it alters the 
critical temperature by a factor of 1 + l/S — 1.4. It is 
even less important in the common circumstance where 
the most important spin fluctuations are collective. This 
classical approximation is, therefore, a less serious source 
of uncertainty in our representation of (III,Mn)V ferro- 
magnet thermodynamics than the modelling issues dis- 
cussed above; in our judgement the substantially greater 
complexity of quantum calculations is not warranted at 
present. Thus, the state of each ion spin Si is specified 
by two polar angles §i, ifi. In interpreting results, how- 
ever, it is important to realize that some low-temperature 
power laws will be altered by quantum freeze out. 

We are interested in thermal expectation values of the 
form 

/ = ~ jf ff ^y"^smtfTr{/>, ¥ 0e-' 3W } , (6) 

where /3 is the inverse temperature, Z the partition func- 
tion, and i?, ip are shorthand notations for the whole set 
of classical spin coordinates. The quantity /($, tp) is a 
function of the ion spin angles and an operator with re- 
spect to the quantum mechanical carrier degrees of free- 
dom over which the trace is performed. In practice we 
replace the fermion trace by a ground state expectation 
value, since the temperatures of interest will always be 
much smaller than the Fermi energy. For typical carrier 
densities n of order O.lnm -3 , the Fermi temperature for 
the carriers is > 1000K, compared to ferromagnetic crit- 
ical temperatures ~ 100K. Therefore, the average over 
the carrier system may be approximated by an average 
over its T = ground state. Thus, 

/ = ! f 7 'dep r cMsm$(0\f(V,^\0)e-^ H ^ , (7) 
^ Jo Jo 

where |0) denotes the ground state of non- interacting 



fermions with the appropriate band Hamiltonian and a 
Zeeman-coupling term h whose effective magnetic field 
B c g is due to exchange interactions with the localized 
spins, 

h = - J d 3 r s(r) ■ B eS (f) 
Bcx(r) = ~Y, J (^- 8i)S&i (8) 

where Cli — (sin(#r) cos(0i),sin(#r) sin(0i), cos(#r)) is 
the direction of the classical spin at Ri. In the follow- 
ing we denote thermal expectation values of quantities 
defined in terms of classical spin orientation variables by 
(•) and quantum mechanical expectation values within 
the carriers ground state by (0| • |0). The latter quanti- 
ties will be averaged also thermally according to Eq. (0). 

III. THE HYBRID MONTE CARLO ALGORITHM 

A standard way to evaluate expectation values of the 
form Eq. (0) is to use classical Monte Carlo algorithms 
which perform a random walk in the phase space of the 
classical variables ($,y>). The probabilities governing this 
Monte Carlo dynamics are specified by the dependence of 
many-fermion energy on the localized-spin configuration. 
The many-fermion ground state is a Slater determinant 
whose single-particle orbitals are the lowest energy eigen- 
states of a single-band or multi-band Hamiltonian of the 
form (|l|). For the case of a parabolic band, the matrix 
elements of the corresponding one-particle Hamiltonian 
in a plane-wave basis read 

(kW\H\ka) = 0%,<W< 

+ £ J ^ ~ k')e^-^'hi ■ , (9) 
i 

where k and a denote wavevector and spin indices, re- 
spectively, J(k) is the Fourier transform of J(r), and L 
the edge length of the simulation cube. Periodic bound- 
ary conditions restrict the admissible values of wavevec- 
tor components to integer multiples of 2tt/L. In Eq. (||) 
t is the vector of Pauli spin matrices. 

Since the many-particle groundstate of the carrier sys- 
tem has to be computed at each Monte Carlo step, the 
computational effort required for the present calculations 
is much larger than in simple classical spin models. In 
the usual Metropolis algorithm, a single spin orientation 
is altered at each step. If this algorithm were employed 
here the time required to diagonalize the single-particle 
Hamiltonian each time would severely limit the efficiency 
of the algorithm. We therefore use the Hybrid Monte 
Carlo algorithm, which was introduced in the mid 1980's 
in the context of lattice field theories . In this method 
all classical variables are altered in one Monte Carlo step. 
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This drastically reduces the number of matrix diagonal- 
isations required to explore statistically important mag- 
netic configurations. The Hybrid algorithm is a powerful 
method for Monte Carlo simulations in systems contain- 
ing coupled classical and quantum mechanical degrees of 
freedom. Very recently, a variant of this algorithm has 
been applied to a system of classical degrees of freedom 
coupled to non- interacting lattice fermions p3[ , a prob- 
lem similar to the one studied here. 

The Hybrid Monte Carlo algorithm determines average 
values defined in terms of a probability distribution P of 
the form 



P{(j>) tx e~ sw 



(10) 



The "action" S depends on a set of classical variables 
summarized by the symbol <f>. The basic trick of the al- 
gorithm is to introduce a "fake dynamics" for the classical 
variables which is governed by the "fake Hamiltonian" 



H' 



1 



(11) 



with "fake momenta" it and an action S' which is not 
necessarily identical to S. 

In this algorithm one Monte Carlo step is performed 
in the following way: (i) Choose a value for each fake 
momentum from its Gaussian distribution, (ii) Let the 
system evolve in Monte Carlo time r according to the 
Hamilton equations of motion, 



d T <j) = 



dW 
dir 



dH' dS' 

C? T 7T = JT— = 7T— ) (12) 



which leads to a new configuration of fields and fake mo- 
menta (</>', 7r'). (iii) Accept this new configuration with a 
probability 



mm{l,exp[ff(<M-#(4>V)]}, 
where H is the "true Hamiltonian" , 



H=-7T 2 



■S(4>). 



(13) 



(14) 



5 = /3(0|W|0)-^ln|sini? 7 



(16) 



The first term is the energy of the carrier system for a 
given classical spin configuration and the logarithms arise 
from the Jacobi determinant in Eq. (Q). Note that this 
part is not multiplied by the inverse temperature (3. 

As pointed out above, the action S 1 need not to be 
the same as S. In our calculation, we use this freedom to 
reduce the computational effort. Our fake action S' is de- 
fined as in Eq. ( |l6| ) but for an itinerant-carrier state |0)j 
that is evaluated in the beginning of each Monte Carlo 
step and taken to be fixed in the susequent integration 
of the fake dynamics (|l^). In addition, we replace J{f) 
by JpdS{r), i.e., we use the ao — > limit of the exchange 
coupling for the fake dynamics. The force term in the 
effective dynamics is then readily evaluated from 



dU(o\n\o)i 
an i 



JS 



(,(01^)10), 



(17) 



In the fake dynamics the localized spins move under the 
influence of the fixed effective magnetic fields generated 
by exchange interactions with the fixed electronic spins of 
|0) j. With this approach, only one matrix diagonalization 
is required to update all Mn spin orientations. However, 
since the acceptance probability (13) is determined by the 
action S, the sampling is still performed with respect to 
the original distribution (|l^) . The choice of S' being dif- 
ferent from S is, therefore, not an approximation as long 
as the Hamiltonian fake dynamics is reversible, a feature 
which can be incorporated in the numerical integration 
using the leapfrog algorithm Jl4] ]. 

The action ( |l6| ) becomes singular if one of the polar 
angles $ is equal to an integer multiple of ir. This may 
lead, in principle, to difficulties in the numerical imple- 
mentation of our method. But since this singularity is 
only logarithmic, and since in the fake dynamics (|l2|) the 
$'s are repelled from the singular points, we have never 
encountered any practical difficulty due to these coordi- 
nate singularities. 



This acceptance condition is reminiscent of the Metropo- 
lis algorithm. We note that the acceptance probability 
is strictly unity for S — S' , and that the acceptance rate 
in this Monte Carlo procedure reflects the difference be- 
tween the two actions S and S' . 

Thermodynamic averages are calculated using 



1 



n=l 



(15) 



where n labels the accepted configurations. As shown in 
Ref. [jl4| this Hybrid Monte Carlo algorithm generates a 
Markov chain which converges to the distribution (|l^) . 

Our calculations adopt this scheme for the set of clas- 
sical spin angles and the action 



IV. NUMERICAL RESULTS 

In this section we present our numerical results. We 
concentrate on the spin polarizations of the manganese 
ions and the carriers as a function of temperature and 
address the ferromagnetic transition. We start with the 
case of parabolic bands described by the Hamiltonian (jl]) . 



A. Parabolic bands 

In the following we present results in dimcnsion- 
ful units for parameters in the range of interest for 
(Ga,Mn)As. It is hoped that ferromagnetism will be re- 
alized in other materials described by the same model 
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but with different parameter values. In this connec- 
tion, it is useful to observe that band and exchange 
interaction terms in our model Hamiltonian both scale 
simply with changes in length scale and model param- 
eters. Results for the parabolic band model have a 
non-trivial dependence only on the ratio of the density 
n band carriers and the concentration of Mn ions N, 
and on the ratio of exchange interaction and band en- 
ergy scales, A/ep- Here A = J p dNS is the continuum 
model mean-field-theory band exchange splitting, and 
e F /k B = Tp » 4230K(m/m*)(n[nm- 3 ]) 2 / 3 is the band 
Fermi energy in the paramagnetic state. There is also 
a very weak dependence on the dimensionless range of 
the exchange interaction clqN 1 / 3 , that does not play an 
important role and will usually be ignored. 

In the following subsections we first show some typical 
numerical results for magnetization averages and fluctu- 
ations and explain how we extract T c from them. We 
then proceed with some technical considerations that en- 
ter into these calculations before summarizing the nu- 
merical results we have obtainted to date. 



1. Magnetisation curves and fluctuations 

Figure |l| shows typical magnetisation data as a func- 
tion of the temperature. These results were obtained 
for a maganese ion density of N = l.Onm" 3 , a carrier 
density n = 0.1nm~ 3 in a cubic simulation volume of 
V = 540nm 3 , i. e. the system contains 540 maganese 
ions and 54 carriers. The effective band mass is half 
the bare electron mass, and the exchange parameter is 
Jpd = 0.15eVnm 3 . The main panel shows the average 
polarisation of the manganese spins, 

M = ^y(\S tot \) , (18) 

i. e. the thermally averaged modulus of the total ion spin, 
along with the carrier magnetisation, 

™ = ^<|<0|.s tot |0)|), (19) 

which is the ensemble average of the modulus of the total 
ground-state carrier spin. Both quantities are divided by 
the number of particles and are close to their maximum 
values at low temperatures. At higher temperatures they 
show the expected transition to a paramagnetic phase. 
The critical temperature of this ferromagnetic transition 
is most readily estimated from numerical results for the 
magnetisation fluctuations: 

3Afn = ^ ? (<|5 ? tot | 2 )-<|5 tot |) 2 ) , (20) 

9p = ^7 ((K°l^l°)| 2 ) - (l<0|s tot |0)|) 2 ) . (21) 

These two fluctuations per particle are plotted in the 
insets of Fig. 111. They both show a pronounced peak 



at a temperature T ~ 100K, defining the finite-system 
transition temperature for these model parameter values. 
In fact, in a region around this transition both datasets 
differ just by a factor of approximately 25, which is the 
square of the ratio of the two spin lengths entering the 
expressions ( p0[) and (^l|), respectively. This observation 
shows explicitly that the correlation length both for the 
manganese ions and in the carrier system is the same near 
the transition, namely given essentially by the system 
size. 

In conclusion, our Monte Carlo approach clearly repro- 
duces the expected ferromagnetic transition. The tran- 
sition temperature T c can be determined unambiguously 
and consistently from the positions of very pronounced 
peaks in total magnetisation fluctuations of both the Mn 
ions and the carriers. 



2. Technical Considerations 

(i) Monte Carlo parameters and disorder real- 
isations. The data of Fig. Oj was obtained by averag- 
ing over five different realisations of the manganese-ions. 
For each Monte Carlo run the system was thermalized 
within the first 1000 accepted steps, and the subsequent 
measurement phase included 10000 accepted steps. The 
Hamilton equations of motions ( O ) are integrated over 
an interval of length 1 in each Monte Carlo step. These 
simulation parameters were used for all results reported 
in this paper. After the first 1000 steps the magnetisa- 
tions for the five different systems differ only by a few 
percent. This indicates that the thermalization phase 
was long enough. In fact, for not too large density ratios 
n/N, the final results for the different disorder realisa- 
tions differ only very weakly from each other. This is 
illustrated in Fig. |^, where the magnetisation curves un- 
derlying the averaged results of Fig. |l| are plotted. Those 
five datasets are hardly distinguishable from each other. 

Finally we mention that possible errors bars due to 
the statistical uncertainty in the magnetisation data for 
individual disorder realisations are by far smaller than 
the symbols used in the diagrams. This follows from 
the analysis of the Monte Carlo magnetisation data as a 
function of Monte Carlo time which shows that the data 
is very well- converged for the algorithm parameters given 
above, an observation which is of course consistent with 
(and actually a necessary condition for) the very smooth 
form of of our data plots as a function of temperature. 

(ii) Wavevector cutoff. Periodic boundary condi- 
tions restrict wavevector components for the carriers to 
integer multiples of 2tt/L, where L = y 1 / 3 is the edge 
length of the simulation cube. In our numerical calcu- 
lations we include wavevectors k up to a certain cutoff 
k c . For the single-parabolic band model the wavevec- 
tor cutoff corresponds simply to a kinetic-energy cutoff. 
The wavevector cutoff has to be chosen carefully, since 
our results are most cleanly interpreted when their cut- 
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off dependence has been saturated. On the other hand 
the dimension of the single-particle Hamiltonian matri- 
ces, whose construction limits our numerical calculations, 
grows with the third power of the cutoff. Fortunately, our 
data converge comparatively quickly with the cutoff pa- 
rameter. Figure [? shows magnetization curves for the 
same system parameters as in Fig. [j] (but in a smaller 
simulation cell of volume V — 140nm 3 ) for two differ- 
ent cutoff parameters. In the first case, k c — 2(2ir/L) 
and the single-particle Hilbert space for the 14 carriers 
has dimension 66 (including the spin degree of freedom) , 
while in the second case, k c = \Z6(2n/L) and the dimen- 
sion increases to 162. The datasets shown in Fig. |3| are 
very close to each other over a large temperature range 
around the ferromagnetic transition. This demonstrates 
that the smaller cutoff is already sufficient (nevertheless 
we used the larger one for Fig. |l|). In the simulations of 
larger systems with more carriers the cutoff parameter is 
appropriately adjusted to a value of k c — y/6(2ir / L) . 

(iii) Finite-size effects. The finite size of the simula- 
tion cell implies that the set of fermion energies obtained 
for a particular Mn spin orientation configuration is dis- 
crete. The gaps between adjacent energy eigenvalues de- 
pends on the size of the simulation cube and may lead to 
systematic fc-space shell effects in the model's finite-size 
dependence. These are less important at temperatures 
of interest, since the typical orientation configuration is 
complex and produces an exchange field that removes all 
degeneracies and smears the fc-space shells. Shell effects 
are therefore more important for weaker exchange cou- 
pling, limiting our ability to do reliable calculations in 
the weak-coupling limit. As we discuss later, however, 
we believe that mean-field theory tends to be reliable in 
this limit, reducing the motivation for Monte Carlo stud- 
ies. Nevertheless, finite-size effects can be observed at all 
coupling strengths. Their character is easily understood 
by recognizing that Mn spins are coupled by the polar- 
ization they produce in the electron system. When a k- 
space shell is partially filled, its degeneracy is lifted by ex- 
change copuling with a Mn spin-configuration. Occupy- 
ing the lower energy states in the split multiplet produces 
a larger electron spin polarization than when the multi- 
plet is full. Of course these effects become less impor- 
tant at larger system sizes when the typical level shift is 
much larger than the typical level spacing. The ferromag- 
netic transition temperature has systematic minima for 
the closed-shell 'magic' electron numbers at which most 
of our simulations are performed. This effect becomes 
smaller with increasing system size. Figure |] shows mag- 
netisation curves for the same densities and system pa- 
rameters as in Fig. [I] for systems with 14, 38, 54, and 30 
carriers. The first three carrier numbers correspond to 
"closed shell" configurations in the paramagnetic phase 
while the last one lies in between. The comparison of 
the data for the three largest carrier numbers shows that 
the finite-size effect has already become weak for carrier 
numbers nV > 30, while the magnetisation curve for the 



system of 14 carriers shows slightly larger values at given 
temperature. However, the finite-size critical tempera- 
ture obtained in the case lies only by about 15K above 
the values for the other three system sizes. 

(iv) The regularisation parameter ciq. As men- 
tioned above we have chosen the regularisation parame- 
ter in the exchange coupling (||) to be ao = O.lnm. This 
seems to be a reasonable value compared with the lat- 
tice constant of GaAs. Figure |^ shows the magnetisation 
data for the same system parameters as in Fig. [j] and a 
simulation volume of V = 140nm 3 for different values of 
ao- The data for ao = O.lnm and ao = 0.2nm are almost 
identical. For larger values of ao a slight departure sets 
in. As a consequence, our results are practically inde- 
pendent of ao over the (physically motivated) range of 
considered values. 

In summary, we have demonstrated that our Monte 
Carlo results are, in a large vicinity of T Cl stable with 
respect to the influence of the wavevector cutoff and the 
finite system size. Moreover, the effect of different disor- 
der realisations and the precise value of the regularisation 
parameter of the exchange coupling is found to be very 
weak. 



3. Results for T c 

We now turn to the transition temperature T c for the 
parabolic band model (|l|). Within mean field theory 
this quantity is given by 



T, 



MF 



Xp 



S(S 



1)NJ? 



pd 



(9*Vb/2) 
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(22) 



where S = 5/2 is the spin length of the manganese ions, 
and g* is the g- factor of the carriers. Their Pauli suscep- 
tibility xp reads 



XP = ~ — (SW2) 2 



(23) 



and is via the Fermi energy e p proportional to the effec- 
tive mass to*. Our objective here is not to make a quan- 
titative prediction of the critical temperature in partic- 
ular ferromagnetic semiconductor systems. In our view, 
uncertainly in model parameters and the possible impor- 
tance of elaborations to the current model that account 
for example for direct Mn ion interactions, make such 
an effort premature. By doing a numerically exact cal- 
culation for a model that captures much of the physics, 
however, we hope to shed light on the range of valid- 
ity and the sense and magnitude of likely corrections to 
mean-field-theory estimates. 

The above expression for T c can be obtained by averag- 
ing the ion-spin and carrier polarization over space. The 
effective field which each manganese spin experiences due 
to a finite carrier polarization is constant in space and 
the carrier bands are in turn spin split by A = J p( iNS 



7 



by a spatially homogeneous effective magnetic field. (The 
limit in which mean-field theory is exact can be achieved 
in our model by letting ao — > oo in Eq. (||).) In this very 
simplified approach, spatial fluctuations and correlations 
between carriers and manganese spins are neglected. As 
a result, mean-field theory predicts T c to be quadratic 
in the exchange parameter J p d and linear in the effective 
band mass m* . This mean- field theory has been used to 
predict values for T c in various semiconductor host ma- 
terials [||. It is, however, expected (l^Jl^l that due to 
the neglect of correlations, Weiss mean-field theory can 
substantially overestimate the critical temperature. 

We therefore investigate the critical temperature T c 
with the help of our Monte Carlo scheme. As seen in 
Fig. ^ this quantity can be read off unambiguously from 
the magnetic fluctuations of both the Mn ions and the 
carriers. To limit computational expense we concentrate 
on systems with fourteen carriers. From Fig. |] we con- 
clude that finite-size effects will lead to values of T c which 
are only slightly too large (less than 20K). Furthermore, 
the qualitative behavior of the T c data as a function of 
system parameters such as the effective mass m* will not 
be affected by the finite simulation-cell size. 

In Fig. |^ we show results for a manganese densities 
N = l.Onm -3 and a mean- field band splitting A = 
JpdNS = 0.5eV. The left panel shows the dependence 
of the critical temperature on the carrier effective mass. 
This dependence is very important for the search for di- 
luted magnetic semiconductor systems with T c 's larger 
than room temperature. In the mean-field approxima- 
tion, T c grows linear with increasing mass. The Monte 
Carlo results clearly deviate from this prediction suggest- 
ing a saturation of T c at carrier masses close to the bare 
electron mass. For even higher masses, we expect the 
electrons to behave more classically and localize around 
individual Mn spins. The cost in electronic energy of 
changing relative orientations of nearby Mn spins will 
get smaller causing T c to decline and eventually causing 
ferromagnetism to disappear. This physics also occurs in 
a continuum model where the s pin stiffness declines as 
1/m* for large band masses [|l2| , |l3| , although there are 
differences in detail. 

To discuss the critical temperature as a function of 
the exchange coupling parameter J p d we observe that the 
Hamiltonian (]l|) satisfies the scaling relation 

m{m*,J v d) = -H\—,qJ p A (24) 

with q > 0. Therefore the saturation of T c as a function 
of the effective mass at fixed J p d corresponds to a linear 
dependence of T c on J p d at fixed m* . This contrasts with 
the mean-field prediction T^f F oc J^ d . 

In the right panel of Fig. we show T c as a function of 
the carrier density. Here the Monte Carlo approach also 
clearly yields a lower critical temperatures than mean- 
field theory. For still higher carrier densities the typical 
distance between nearby Mn ions will become larger than 



the band electron Fermi wavelength, causing the sign of 
the typical exchange coupling to oscillate in an RKKY 
fashion. We expect the resulting frustration to make the 
ferromagnetic state unstable, possibly leading to a regime 
of spin-glass order. 

As discussed above, different disorder realisations with 
respect to the ion positions give almost identical critical 
temperatures (see also Fig. |2|). This does not mean, how- 
ever, that the presence of disorder is unimportant. Nu- 
merical finite-size studies at T = have shown that in 
the limit of strong exchange couplin g J pt j the presence of 
disorder increases the spin stiffness |1^J2J] for the regime 
of densities and masses covered by our calculations. This, 
in turn, enhances the critical temperature since, accord- 
ing to analytical estimates from spin- wave theory, T c is 
approximately proportional to the spin stiffness uM . On 
the other hand, our Monte Carlo results show that the 
spin stiffness enhancement does not depend on the con- 
crete disorder realization. 



4- Local vs. global polarisation 

In mean-field theory, the free carrier band spin- 
splitting vanishes as the critical temperature is ap- 
proached. This is, however, not generally the correct 
physical picture. When long-wavelength collective fluctu- 
ations Jl^,[l3| , which are neglected by mean-field theory, 
drive the breakdown of ferromagnetism the local carrier 
spin polarization can remain finite even above T c . Fer- 
romagnetism and the global spin polarization disappears 
only because of the loss of long-range spatial coherence. 
The applicability of this picture to the studied parameter 
ranges of the kinetic exchange model is confirmed by our 
Monte Carlo studies. In Fig. we compare the global 
(Eq. (|i~9|)) and local carrier spin polarizations, 

which is the thermally and spatially averaged ratio of 
the modulus of the carrier spin density s(r ) and the local 
carrier density n(f ). The system parameter are the same 
as m Fig. 0. While the global spin polarization vanishes 
at the critical temperature, the local polarization remains 
finite and saturates at about fourty percent of its T = 
value. 



B. Six- band model 

We now turn to the six-band model in which the kinetic 
energy part is given by the Kohn-Luttinger Hamiltonian 
(ff). To model GaAs we use (71,72,73) = (6.85,2.1,2.9) 
for the Luttinger parameters, and the spin-orbit coupling 
energy is A so = 0.34eV. 

Figure || shows typical magnetisation data for a system 
with exchange coupling J p d — 0.15eVnm 3 , carrier density 
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is n = 0.1nm~ 3 , and Mn ion density N = 1.0nm~ 3 in a 
volume of V — 280nm 3 . The wavevector cutoff here is 
k c = 2(2n/L) offering 192 single-particle states to the 28 
carrier in the system. 

As in the case of parabolic bands, a ferromagnetic tran- 
sition is clearly signalled by pronounced peaks in the 
magnetic fluctuations of both Mn ions and carriers. We 
find that, in contrast to the parabolic two-band model, 
the carrier magnetisation is already reduced at tempera- 
tures well below T c . Another difference compared to the 
parabolic-band model concerns the shape of the magnetic 
fluctuations for the manganese ions and the carriers as a 
function of temperature. Although both curves indicate 
the same value for T c , their shape in the vicinity of T c is 
slightly different and the ratio of the their fluctuations is 
smaller than 25 (the square of the ratio of spin lengths 
involved). These differences arise because of the more 
complicated band structure and the spin-orbit coupling 
present in the Kohn-Luttinger Hamiltonian. 

As for the two-band model, the data shown in Fig. ||are 
the average over five different disorder realisations. Also 
in the present case, the dependence of the data on the 
concrete realisation is very weak. Moreover, we again find 
that the local carrier polarization remains finite above T c 
while the global magnetization vanishes. 

In the right panel Fig. || we plot the transition temper- 
ature as a function of the exchange coupling J pc i for two 
different system sizes. Both data sets agree within error 
bars and show a linear dependence of T c on J pc i- This 
finding is the same as for the two-band model and con- 
trasts with mean- field theory which predicts T c cx J^ d . 
The left panel of Fig. ^ shows the transition temperature 
as a function of J pc i for the parabolic model for an effec- 
tive mass m* — 0.5m. This value is close to the heavy- 
hole mass in the Kohn-Luttinger model for parameters 
appropriate for GaAs as given above. The data in the 
left panel can be obtained from the left panel of Fig. || 
via the scaling relation (p4|). Comparing the two pan- 
els of Fig. ^ demonstrates that, in the range of carrier 
densities studied, a single parabolic band with an effec- 
tive mass close to that of the heavy-hole Kohn-Luttinger- 
model band provides a reasonably good approximation to 
the behavior of the six-band system. 

V. SUMMARY 

We have given a detailed report on Monte Carlo stud- 
ies of a model that captures essential aspects of ferro- 
magnetism in (III,Mn)V ferromagnets. In the kinetic- 
exchange model we study localized magnetic moments 
formed by Mn ions are coupled antiferromagnetically to 
carriers (holes) in the valence band. In most of our sim- 
ulations we describe the valence-band carriers by simple 
parabolic bands but we have investigated a more realistic 
six-band k-p model as well. We treat the manganese spins 
as classical magnetic moments, and take the carriers ef- 



fectively at zero temperature. Both approximations are 
comparatively weak. Our treatment accounts for spa- 
tial fluctuations in the magnetization and for disorder 
induced by the randomness of the Mn ion positions, ef- 
fects which are both effects are neglected in a mean-field 
description. 

We outline the Hybrid Monte Carlo algorithm em- 
ployed and discussed its advantages compared to other 
methods. A complete account of important technical de- 
tails of the computation is given. We have shown that 
the effects of different disorder realizations, wavevector 
cutoff, and finite-size effects are well under control. Our 
results depend weakly on the regularization parameter clq 
which describes the spatial range of the exchange mech- 
anism between carriers and ions. 

We find that due to spatial fluctuations and correla- 
tions the total magnetization and therefore also the fer- 
romagnetic transition temperature T c is considerably re- 
duced in comparison to mean-field estimates, increasingly 
so as the exchange coupling strength and the band ef- 
fective mass are increased. In contrast with mean-ficld- 
theory predictions, T c is not a linear function of the ef- 
fective band mass and is not proportional to the square 
of the exchange coupling constant. The discrepancy be- 
comes more severe for higher predicted mean-field critical 
temperatures. It remains an interesting to what extend 
these deviations from mean field theory arc influenced 
by the presence and type of disorder with respect to the 
manganese positions. 

The importance of spatial fluctuations of the magneti- 
zation, neglected in a mean-field description, is empha- 
sized by the fact that the local carrier spin polarization 
remains finite even well above the ferromagnetic transi- 
tion temperature. 
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FIG. 1. Magnetization curves for manganese ions and carri- 
ers. The upper and lower inset show the magnetic fluctuations 
for the manganese ions and the carriers, respectively. Both 
differ by a factor of approximately 25 reflecting the square of 
the ratio of spin lengths. The density of manganese ions is 
N = l.Onm" 3 , the carrier density is n = O.lnm 3 in a cubic 
volume of V = 540nm 3 . The band mass is half the bare elec- 
tron mass with an exchange parameter of J p d = 0.15eVnm 3 . 
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FIG. 2. Magnetisation curves for five different realisations 
of manganese positions underlying the averaged data of Fig. 
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FIG. 3. Magnetisation data for the same system parame- 
ters as in Fig. hi in a cubic volume of V = 140nm 3 for two dif- 
ferent values of the wavevector cutoff. In a large area around 
the ferromagnetic transition both datasets are very close two 
each other. 
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FIG. 5. Magnetisation data for the same system parame- 
ters as in Fig. hi in a cubic volume of V = 140nm 3 for different 
values of the regularisation parameter ao. At not too large 
values for ao the dependence on this quantity is extremely 
weak. 
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FIG. 4. Magnetisation data for the same densities and 
band mass and exchange parameter as in Fig. [l] for differ- 
ent system sizes. The first three systems with 14, 38, and 
54 carriers correspond to "closed shell" configurations in the 
paramagnetic carrier state, while the last system (30 carriers) 
lies in between. 
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FIG. 6. The critical temperature T c as a function of the 
carrier mass (left panel) and the carrier density (right panel) . 
The exchange parameter is in both cases J p d = 0.20eVnm 3 
leading to a zero temperature mean-field spin splitting of 
A = J pd NS = 0.5eV. The results of the Monte Carlo runs 
are compared with the mean-field predictions. 
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FIG. 7. Local versus global carrier spin polarisation for the 
same system as in Fig. hi as a function of temperature. The 
local carrier polarisation remains finite above T c . 
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FIG. 8. Magnetisation curves for Mn ions and carri- 
ers in the six-band model for n exchange coupling of 
J pd = 0.15eVnm 3 . The carrier density is n = 0.1 nm 
with an Mn ion density of TV = l.Onm -3 in a volume of 
V = 280nm 3 . As in the case of parabolic bands, the fer- 
romagnetic transition in clearly and consistently signalled by 
pronounced peaks in the magnetic fluctuations shown in the 
insets. 
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FIG. 9. The right panel shows critical temperature T c as a 
function of the exchange parameter J p d for the same particle 
densities as in Fig. ^|and two different system sizes. Both data 
sets agree within error bars and show a linear dependence 
of T c on Jpd. In the left panel, the corresponding data for 
the parabolic two-band model with an effective mass of half 
the bare electron mass is plotted. The latter system is a 
reasonable approximation to the six-band case. 
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